#########################################
# RUN THE IV MODELS WITH ALTERNATIVE DV #
#########################################

# Author: Kasia Nalewajko
# First created: 13 December 2023
# Replicated: 12 June 2024

rm(list = ls())

# LOAD PACKAGES -----------------------------------------------------------

if (!require("dplyr")) install.packages("dplyr")
if (!require("estimatr")) install.packages("estimatr")
if (!require("modelsummary")) install.packages("modelsummary")
if (!require("fixest")) install.packages("fixest")
if (!require("sf")) install.packages("sf")
if (!require("spdep")) install.packages("spdep")

# LOAD DATA ---------------------------------------------------------------

load("./00 SUBMITTED/00 APSR final/04 replication_files/01 data/main.Rda")

# RUN MODELS --------------------------------------------------------------

# Model with limited controls (columns 1 and 2) ----

mlim_county <- fixest::feols(log(all_jewish_victims+1) ~ log(pop1936+1) + synagogues_dummy + log(collabos_antijew1000+1) + log(DHI_milipol_sum1942+1) + area_sqkm  + longitude + longitudesq + latitude + latitudesq | zone5 + as.factor(ar_name)
                             | log((FFIFFCgentile*1000/pop1936)+1) ~ log(mpf1000+1),
                             cluster = ~c(id_bureau, resregion2),
                             data = main)

# extract additional statistics
finfo <- fitstat( mlim_county, "ivf")
mlim_county_fstat <- round(finfo[["ivf1::log((FFIFFCgentile * 1000/pop1936) + 1)"]][["stat"]])
mlim_county_wh <- round(summary(mlim_county)$iv_wh[["p"]], digits = 4)

# run Moran test and extract the statistic
cantons_sf <- sf::st_read("./data/other/maps/france1940/CANTONS_1940/CANTONS_1940.shp")
obsRemoved <- abs( mlim_county[["obs_selection"]][["obsRemoved"]])
cantons_sf <- cantons_sf[-obsRemoved,]
spatial_weights_matrix <- poly2nb(pl = cantons_sf, row.names = cantons_sf$deppct, queen = T)
spatial_weights_matrix <- nb2listw(spatial_weights_matrix, zero.policy = T) 
mlim_county_res <- residuals( mlim_county)
mlim_county_moran <- moran.test(x = mlim_county_res, listw = spatial_weights_matrix, zero.policy = T)
mlim_county_moran <- round( mlim_county_moran$estimate[[1]], d = 3)


# Fully specified model (columns 3 and 4) ----

mcontrols_county <- fixest::feols(log(all_jewish_victims+1) ~ log(pop1936+1) + synagogues_dummy + catholic_churches + log(collabos_antijew1000+1) + log(DHI_milipol_sum1942+1) + log(FRANCISTE_36_1_pct+1) + log(ActionFrancaise_19_pct+1) + turnout_36_1_pct + log(nuance_centredroit_36_1_pct+1) + log(nuance_droite_36_1_pct+1) + log(nuance_centregauche_36_1_pct+1) + log(nuance_gauche_36_1_pct+1) + log(nuance_extremegauche_36_1_pct+1) + area_sqkm  + longitude + longitudesq + latitude + latitudesq | zone5 + as.factor(ar_name) #  + log(nuance_extremedroite_36_1_pct+1)
                                  | log((FFIFFCgentile*1000/pop1936)+1) ~ log(mpf1000+1),#  + turnout_14_pct 
                                  cluster = ~c(id_bureau, resregion2), 
                                  data = main)

# extract additional statistics
finfo <- fitstat( mcontrols_county, "ivf")
mcontrols_county_fstat <- round(finfo[["ivf1::log((FFIFFCgentile * 1000/pop1936) + 1)"]][["stat"]])
mcontrols_county_wh <- round(summary( mcontrols_county)$iv_wh[["p"]], digits = 4)

# run Moran test and extract the statistic
cantons_sf <- sf::st_read("./data/other/maps/france1940/CANTONS_1940/CANTONS_1940.shp")
obsRemoved <- abs( mcontrols_county[["obs_selection"]][["obsRemoved"]])
cantons_sf <- cantons_sf[-obsRemoved,]
spatial_weights_matrix <- poly2nb(pl = cantons_sf, row.names = cantons_sf$deppct, queen = T)
spatial_weights_matrix <- nb2listw(spatial_weights_matrix, zero.policy = T) 
mcontrols_county_res <- residuals( mcontrols_county)
mcontrols_county_moran <- moran.test(x = mcontrols_county_res, listw = spatial_weights_matrix, zero.policy = T)
mcontrols_county_moran <- round( mcontrols_county_moran$estimate[[1]], d = 3)

# Preview ----

modelsummary(list(
  'First Stage' = mlim_county$iv_first_stage[[1]],
  'Second Stage' = mlim_county,
  'First Stage' = mcontrols_county$iv_first_stage[[1]],
  'Second Stage' = mcontrols_county),
  stars = c('*' = .1, '**' = .05, '***' = .01))

# OUTPUT TO LATEX ---------------------------------------------------------------

# add_rows
add_stats <- tibble::tribble(
  ~x, ~y, ~q, ~w, ~l,
  'F stat. (1st stage)',  mlim_county_fstat,  mlim_county_fstat,  mcontrols_county_fstat,  mcontrols_county_fstat, 
  'Moran stat.',  mlim_county_moran,  mlim_county_moran,  mcontrols_county_moran,  mcontrols_county_moran,
  'Wu-Hausman p-value',  mlim_county_wh,  mlim_county_wh,  mcontrols_county_wh,  mcontrols_county_wh 
)

attr(add_stats, 'position') <- c(39, 40, 41)

msummary(list(
  'First Stage' = mlim_county$iv_first_stage[[1]],
  'Second Stage' = mlim_county,
  'First Stage' = mcontrols_county$iv_first_stage[[1]],
  'Second Stage' = mcontrols_county
),
stars = c('*' = .1, '**' = .05, '***' = .01),
gof_omit = "AIC|BIC|RMSE|R2 Within|R2 Within Adj.",
add_rows = add_stats,
output = "latex",
coef_rename = c("log(mpf1000 + 1)" = "WWI military death rates",
                "fit_log((FFIFFCgentile * 1000/pop1936) + 1)" = "Insurgent presence",
                "months_Verdun_Petain" = "Months Vérdun Pétain",
                "log(pop1936 + 1)" = "1936 population",
                "log(collabos_antijew1000 + 1)" = "Collaborators",
                "log(DHI_milipol_sum1942 + 1)" = "1942 state presence",
                "synagogues_dummy" = "Synagogues",
                "area_sqkm" = "Area size (km2)",
                "longitude" = "Longitude",
                "latitude" = "Latitude",
                "longitudesq" = "Longitude (sq)",
                "latitudesq" = "Latitude (sq)",
                "log(allocated_jpop + 1)" = "1941 Jewish population",
                "catholic_churches" = "Catholic churches",
                "log(FRANCISTE_36_1_pct + 1)" = "Franciste vote 1936",
                "log(ActionFrancaise_19_pct + 1)" = "Action Française vote 1919",
                "turnout_36_1_pct" = "Turnout 1936",
                "log(nuance_droite_36_1_pct + 1)" = "Right vote 1936",
                "log(nuance_centredroit_36_1_pct + 1)" = "Centre-right vote 1936",
                "log(nuance_centregauche_36_1_pct + 1)" = "Centre-left vote 1936",
                "log(nuance_gauche_36_1_pct + 1)" = "Left vote 1936",
                "log(nuance_extremegauche_36_1_pct + 1)" = "Extreme left vote 1936"
),
title = "Main results. Dependent variable: Number Holocaust victims (logged)"
)

